#!/usr/bin/python
# GM_alterpars.py

import math
import numpy
from time import strftime

import GM_InterpRoutines as IR
import GM_Swap

secperyear = 3600*24*365.25
kmperpc = 3.085e13

def alterpars(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    tiny = 1e-6
    # Fix units, types for simpars
    isacluster = bool(isacluster)
    nodes = numpy.array(nodes, dtype=int)
    scale = float(scale)*1000. # kpc -> pc
    galactocentric = bool(galactocentric)
    Bochanski2010IMFandLF = bool(Bochanski2010IMFandLF)
    BoBand = bool(BoBand)
    system_or_single = str(system_or_single)
    LensedMag = float(LensedMag)
    SourceMag = numpy.array(SourceMag, dtype=float)
    Magnification = float(Magnification)
    brightlimit = int(brightlimit)
    lres = int(lres)
    bres = int(bres)
    dres = int(dres)
    if band == 'V':
        c0 = float(c0_V)
        rho0 = float(rho0_V)
    if band == 'I':
        c0 = float(c0_I)
        rho0 = float(rho0_I)
    fitcrho = bool(fitcrho)
    fixedrho = bool(fixedrho)
    fixedc = bool(fixedc)
    latitudeweight = float(latitudeweight)
    if latitudeweight > 1.:
        latitudeweight = 1.
        print "The value you chose for latitudeweight was too large.  It has been corrected to the maximum value, 1."
    if latitudeweight < 0.:
        latitudeweight = 0.
        print "The value you chose for latitudeweight was too small.  It has been corrected to the minimum value, 0."
    interpolateMv = bool(interpolateMv)
    intervals = float(intervals)
    ApparentLimits = numpy.array(ApparentLimits, dtype=float)
    if ApparentLimits.min() >= SourceMag.min():
        print "Warning!  It is best if the lower boundary of ApparentLimits is less than the lowest value of SourceMag."
    if ApparentLimits.max() <= SourceMag.max():
        print "Warning!  It is best if the upper boundary of ApparentLimits is greater than the highest value of SourceMag."
    interpolateMass = bool(interpolateMass)
    maxmassinterval = float(maxmassinterval)
    timestamp = strftime("%Y%m%d_%H%M%S")
    codeindex = -1
    pencilBeam = bool(pencilBeam)
    pencilBeamPointing = numpy.array(pencilBeamPointing, dtype = float)
    # Note that all the stuff about generating location/coordinate arrays has moved to GM_Locations
    # Fix units, types for phypars
    # Han values
    R0 = float(R0)*1000. # kpc -> pc
    hR = float(hR)*1000. # kpc -> pc
    hdz = float(hdz)
    ScaleList = numpy.array(ScaleList, dtype=float)
    bulk_vels = numpy.array(bulk_vels, dtype=float)*secperyear/kmperpc # km/s -> pc/yr
    sigma_vels = numpy.array(sigma_vels, dtype=float)*secperyear/kmperpc # km/s -> pc/yr
    # Derived parameters, Han, pt.1
    vanillaScaleList = numpy.array(ScaleList)
    # For modifiedScaleList see below
    # Binney & Merrifield values
    MvList = numpy.array(MvList)
    LumFuncList = numpy.array(LumFuncList, dtype=float)
    # Bochanski values
    # First 2 already in correct format
    BoCorrectedValues = dict(BoCorrectedValues)
    for i in BoCorrectedValues:
        BoCorrectedValues[i] = float(BoCorrectedValues[i])
    BoUncorrectedValues = dict(BoUncorrectedValues)
    for i in BoUncorrectedValues:
        BoUncorrectedValues[i] = float(BoUncorrectedValues[i])
    BoSystemLF_rband = numpy.array( numpy.array( BoSystemLF_rband).T)
    BoSingleLF_rband = numpy.array( numpy.array( BoSingleLF_rband).T)
    BoSystemLF_Jband = numpy.array( numpy.array( BoSystemLF_Jband).T)
    BoSingleLF_Jband = numpy.array( numpy.array( BoSingleLF_Jband).T)
    BoSystemMF = numpy.array( numpy.array( BoSystemMF).T)
    BoSingleMF = numpy.array( numpy.array( BoSingleMF).T)

    if Bochanski2010IMFandLF == True:
        if system_or_single == 'single':
            BoMassFunc = BoSingleMF
            if BoBand == 'Bo_r':
                BoLumFunc = BoSingleLF_rband
            if BoBand == 'Bo_J':
                BoLumFunc = BoSingleLF_Jband
        if system_or_single == 'system':
            BoMassFunc = BoSystemMF
            if BoBand == 'Bo_r':
                BoLumFunc = BoSystemLF_rband
            if BoBand == 'Bo_J':
                BoLumFunc = BoSystemLF_Jband
        # Major support needed for Bochanski mass function
        # Additional support for MF/LF choice may be useful
        MvList = BoLumFunc[0,:]
        LumFuncList = BoLumFunc[1,:]*2.e-3
        BoMasses = 10**BoMassFunc[0,:]
        BoMF = BoMassFunc[1,:]*1.e-3
#     = numpy.array()
    # Derived parameters, B & M _OR_ Bochanski
    (delMvList, temp) = IR.DeltasLimitsLinear(MvList, False)
    vanillaMvList = numpy.array(MvList)
    vanillaLumFuncList = numpy.array(LumFuncList)
    vanilladelMvList = numpy.array(delMvList)
    modifiedMvList = IR.NewEvenArray(MvList, intervals**-1, True, False)
    modifiedLumFuncList = IR.OldToNewLinear(MvList, LumFuncList, modifiedMvList)
    (modifieddelMvList, temp) = IR.DeltasLimitsLinear(modifiedMvList, False)
    # Derived parameters, Han, pt.2
    modifiedScaleList = IR.OldToNewLinear(vanillaMvList, vanillaScaleList, modifiedMvList)
    # Cox values
    bvals = numpy.array(bvals)
    MvVals = numpy.array(MvVals)
#    if MvVals.min() > SourceMag.min():
#        print "Warning!  It is best if the lower boundary of MvVals is less than or equal to the lowest value of SourceMag."
#    limitset = numpy.array(limitset)
    if band == 'V':
        N_bM = numpy.array(N_bM_V)
        allsky = numpy.array(allsky_V)
    if band == 'I':
        N_bM = numpy.array(N_bM_I)
        allsky = numpy.array(allsky_I)
    # Derived parameters, Cox
    (delMvVals, limitset) = IR.DeltasLimitsLinear(MvVals, False)
    vanillaMvVals = numpy.array(MvVals)
    vanilladelMvVals = numpy.array(delMvVals)
    vanillalimitset = numpy.array(limitset)
    fittingMvVals = IR.NewEvenArray(MvVals, intervals**-1, True, False)
    (fittingdelMvVals, fittinglimitset) = IR.DeltasLimitsLinear(fittingMvVals, False)
    modifiedMvVals = IR.NewEvenArray(ApparentLimits, intervals**-1, True, True)
    (modifieddelMvVals, modifiedlimitset) = IR.DeltasLimitsLinear(modifiedMvVals, False)

#    print "vanillalimitset is", vanillalimitset
#    print "fittinglimitset is", fittinglimitset

    # Kroupa
    masslimits = numpy.array(masslimits, dtype=float)
    massindices = numpy.array(massindices, dtype=float)
    massmults = numpy.array(massmults, dtype=float)
    MvMass = numpy.array(MvMass, dtype=float)
    # Derived values, Kroupa
    (delMvMass, limitsMvMass) = IR.DeltasLimitsLinear(MvMass, True)
    vanillaMvMass = numpy.array(MvMass)
    vanilladelMvMass = numpy.array(delMvMass)
    vanillalimitsMvMass = numpy.array(limitsMvMass)
    tempMvOnly0 = MvMass[:,0]
    tempMassOnly0 = MvMass[:,1]
    tempMassOnly1 = IR.FillInGaps(tempMassOnly0, maxmassinterval)
    tempMvOnly1 = IR.OldToNewLinear(tempMassOnly0, tempMvOnly0, tempMassOnly1)
    modifiedMvMass = numpy.zeros((len(tempMvOnly1),2))
    modifiedMvMass[:,0] = tempMvOnly1
    modifiedMvMass[:,1] = tempMassOnly1
    (modifieddelMvMass, modifiedlimitsMvMass) = IR.DeltasLimitsLinear(modifiedMvMass, True)
    # Derived values, Han + Kroupa
    MvMassScales = IR.OldToNewLinear(MvList, ScaleList, MvMass[:,0])
    vanillaMvMassScales = numpy.array(MvMassScales)
    modifiedMvMassScales = IR.OldToNewLinear(MvList, ScaleList, modifiedMvMass[:,0])
    # Sofue
    SofueRotV = numpy.array(SofueRotV)
    SofueRotV[:,0] = SofueRotV[:,0]*1000 # kpc -> pc
    SofueRotV[:,1] = SofueRotV[:,1]*secperyear/kmperpc #km/s -> pc/yr
    # Derived values, Sofue
    # I'm not sure if this is something that ought to happen here
#    modifiedSofueRotV = IR.OldToNewLinear( SofueRotV[:,0], SofueRotV[:,1],  "newx")
    # Bring phypars back to inform simpars again
    # NOTE, this usage is fine/expected even in the pencilbeam case
    hzmean = ScaleList.mean()
    ScaleZ = scale*hzmean/hR
    PtScaled = (scale + R0)/dres
    PtScalela = math.pi/lres
    PtScaleba = math.pi*0.5/bres

############

    # An external weight, sinusoidal, that the user can input
    temp = latitudeweight/2.
    latitudeweightexternal = 1 - temp - temp*numpy.cos(bvals*math.pi/90.)
    # An internal weight, natural to the system, based on the angular areas represented by each point
    (temp1, temp) = IR.DeltasLimitsLinear(bvals, True)
    temp = temp*math.pi/180.
    latitudeweightinternal = numpy.sin(temp[:,1]) - numpy.sin(temp[:,0])
#     = bool()
#     = int()
#     = float()
#     = numpy.array()
#     = numpy.array(, dtype=int)
#     = numpy.array(, dtype=float)
    # Here switch to appropriate (modified) variable sets if called for
    # For this section to work properly, gotta bundle in (and maybe bundle out) phypars
    for i in phypars:
        cmd = "phypars['%s'] = %s" % (i,i)
        exec cmd
    # Note: Some additional swapping must go on within fitc0rho0
    if interpolateMv == True:
        phypars = GM_Swap.swap_mags(phypars, False)
    if interpolateMass == True:
        phypars = GM_Swap.swap_mass(phypars, False)

    # Pack it away and ship it out
    for i in simpars:
        cmd = "simpars['%s'] = %s" % (i,i)
        exec cmd
#    for i in phypars:
#        cmd = "phypars['%s'] = %s" % (i,i)
#        exec cmd
    return (simpars, phypars)
